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Abstract 

We report on numerical solutions of the Gross-Pitaevskii equation for two-dimensional flow of a su- 
pcrfluid condensate through a small orifice. Above a critical velocity of about 30% of the speed of sound, 
cavitation occurs in the throat of the orifice. The cavitating bubbles form the cores of singly quantized 
vortices which detach from the boundary and are convected downstream. 
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1. Introduction 

Boundary layer separation is a familiar phenomenon in classical fluid dynamics [1]. ft is invariably 
associated with the appearance of vorticity in the bulk of the fluid. The main flow and the relatively 
stagnant fluid in the lee of the bounding body may be in irrotational motion but they are separated by 
the detached boundary layer, a sheet of fluid containing vorticity generated by viscous forces while it was 
in contact with the body. For example, when a classical fluid is forced through a small opening the flow 
separates close to the narrowest point of the orifice and forms a jet penetrating into the ambient fluid. This 
jet is bounded by a tube of vorticity which confines the flow much as a solenoid confines the magnetic field it 
produces. When the opening has a sharp salient edge the boundary layer separation takes place at arbitrarily 
low flow velocities. If the edge is rounded there will be some critical velocity below which the flow remains 
attached and potential. 

Similar regimes occur when a superfluid is forced through a small opening. At low velocities the superflow 
is potential and dissipationlcss. Liquid can pass through the orifice without the need for a pressure difference 
between the two sides. Above a critical flow rate a pressure head is needed to maintain the flow. At sufficiently 
high flow rate even a superfluid behaves as a classical liquid, and in this regime we can be confident that 
a collimated jet has formed. Because vorticity is quantized in a superfluid, the flow separation requires the 
shedding of discrete vortices from the boundary. These vortices will approximate the classical vortex sheet 
bounding the the jet, just as the individual windings approximate an ideal solenoid. 

In recent years resonators containing micron sized orifices have been constructed in which one can 
monitor the dissipation produced by the shedding of a single vortex [2,3]. At very low temperatures the 
critical velocities in the orifice become independent of temperature, and this temperature independence is 
cited as evidence that the vortices are being created by quantum mechanical tunneling [4,5]. This inference 
is probably correct, but it is worth appreciating that thermal or quantum nucleation are not the only 
alternatives for creating a vortex. Vortices may also be shed from a boundary via cavitation. An accurate 
estimate of the critical velocity for this "classical" vortex creation mechanism is needed before one can state 
with confidence that one is seeing macroscopic quantum tunneling events. 

This present paper provides a rough estimate of this critical velocity. It reports numerical solutions 
of the time-dependent Gross-Pitaevskii (G-P) equation [6,7] for flow through (two dimensional) orifices. 
These solutions demonstrate that when the velocity near the classical flow separation point exceeds about 
30% of the speed of sound, low pressure in the throat of the orifice causes cavitation. Pairs of bubbles are 
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pulled off the boundary, the pair forming the core of what in three dimensions would be a singly-quantized 
smoke-ring-like vortex. 

It is surprising that this classical cavitation effect was not discovered long ago. It is not described in 
any of the standard references on the G-P equation but, while we were performing the calculations reported 
here, our attention was drawn to the recent work of Frisch, Pomeau and Rica who observed similar vortex 
shedding from the boundary of cylinders immersed in a uniform flow [8,9]. These authors interpret the vortex 
shedding as a substitute for shock-wave solutions (which do not occur in the G-P equation) rather than, 
as seems more natural to us, an emulation of boundary layer separation. Also they do not exhibit detailed 
pictures of the formation of the vortices. Consequently, in addition to readvertizing the phenomenon and 
estimating the critical velocity, we include some detailed plots of how the order parameter phase evolves 
as the vortices form. In the second section we will review the well-known interpretation of the real and 
imaginary parts of the Gross-Pitaevskii equation as the equations of fluid mechanics. We also include a 
discussion of the Joscphson phase winding between the core of a jet and the ambient fluid. In the third 
section we will report on the numerical methods and boundary conditions used to induce the flows. In the 
fourth section we describe the results. 

2) The Gross-Pitaevskii equation 

The Gross-Pitaevski (G-P) non-linear Schrodinger equation [6,7] provides a simple model for the dy- 
namics of the condensate in 4 He. Although the local version of this equation does not describe the roton part 
of the quasiparticle spectrum [9], it does model the most essential ingredients: long-range phase coherence 
coupled with mass and momentum conservation. 

To establish our notation, and to develop some physical insight, we begin by reviewing the close re- 
lationship between the G-P equation and classical fluid dynamics. Because it is important to distinguish 
between quantum and classical effects we will include explicit factors of h. For generality we will also include 
a coupling to a gauge field, although the application in this paper is primarily to neutral superfluid 4 He. 

The Gross-Pitaevskii equation is the non-linear Schrodinger equation of the form 



The sign of A will be understood to be positive, so the last term in (1.1) represents a repulsive force between 
the the 4 He atoms. (This is the opposite sign from that usually taken when studying the one-dimensional 




(1.1) 
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non-linear Shrodinger equation as an example of an integrable soliton system.) The mass m should be taken 
to be that of a helium atom (6.647 x 10~ 27 kg). The values of A, and the number density po are usually 
chosen to fit the measured density (mpo = 145kgm~ 3 ), and speed of sound (c = ^Xpo/m = 230ms -1 ) in 
the fluid. 

By using the Madelung transformation [10], (1.1) can be recast as the equation of motion of a charged 
compressible fluid having equilibrium particle- number density p - We set 4 f = y/pe l9 and define a velocity 
field v in such a way that the number-current, 

j = ^— (**(V - ieA/K)^ - ((V + ieA/h)^*)^) , (1.2) 

A 771 % 

may be written as j = pw. This requires 

v= — (Ve-eA/h). (1.3) 
m 

In the absence of vortex singularities in 'J, the vorticity, u) = V A v, is completely determined by the 
gauge field to be w = — — V A A = — — , i.e 

muj + eB = 0. (1.4) 



For neutral superfluids (1.4) implies irrotational motion and hence, at low velocities, where the effects 
of compressibility can be ignored, leads to D'Alembert's paradox (the absence of drag forces). For charged 
superfluids equation (1.4) is responsible for the Meissner effect: A penetrating uniform B field would require 
uniform vorticity, i.e, rigid rotation. A rigidly rotating body possesses a kinetic energy which grows faster 
than the volume of the system and so is impossible in the thermodynamic limit. Alternatively, taking the 
curl of the Maxwell equation VAB = ej and using j = pv s together with (1.4), implies that 

2 

V 2 B-^B = 0, (1.5) 

m 

which leads to flux screening. Note that h does not appear in (1.4) or in the screening length (e 2 p/m) -1 / 2 . 

With the definition (1.3), the imaginary and real parts of (1.1) become, respectively, the continuity 
equation 

d t p + V • pw = 0, (1.6) 
and the Euler equation governing the flow of a barotropic fluid 

m(d t v + v • Vv) = e(E + v A B) - V/x. (1.7) 
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The word barotropic refers to the simplifying property that the pressure term -VP, which occurs on the 
right hand side of the conventional Euler equation, is here combined into the gradient of a potential 

^-ri-^nr- (L8) 

The potential p, contains the expected compressibility pressure, depending on the deviation from the equilib- 
rium density, together with a correction depending on gradients of p. This correction, the quantum pressure, 
is the only place that h appears in the flow equations. It sets the length scale £ = h/y/2m\p over which 
the superfluid density heals after being forced to zero by a boundary or a vortex singularity. With the 
parameters chosen as suggested above we find £ = .487 A. (h also manifests itself, of course, in the quantum 
of circulation, k = 2Trh/m. ) 

The Euler equation (1.7) is derived by first taking the gradient of (1.1) and interpreting the result as 
the Bernoulli equation, 

m(d t v - v A u) = e(E + v A B) - V Qmv 2 + p)j , (1.9) 

which is equivalent to (1.7). 

In (1.9) a cancellation of the mv A u term against the ev A B term is evident upon use of (1.4). It 
is after this cancellation, and so without reference to either B or ui, that the hydrodynamic picture of 
superconductivity is conventionally displayed. It seems preferable to keep both u> and B in (1.9) and rewrite 
it as (1.7). By doing this one can see that the only difference between superfluid dynamics and classical fluid 
dynamics lies in the constraint (1.4). 

In the next section we solve (1.1) numerically for the case of a neutral superfluid forced to flow through 
an orifice. Let us first consider a crude model of the extreme case where the downstream flow forms a tubular 
jet penetrating into superfluid at rest. 

We suppose that there is a pressure head, and correspondingly a difference pi — p m the potential, 
between the asymptotic parts of the reservoirs communicating via the orifice. This means that the phase of 
the order parameter of the fluid at rest in the upstream reservoir, 9\, must be falling behind the phase of 
the order parameter in of the fluid at rest in the lower presure container, 9q, at a rate given by 

-n^(0i - e ) = m - tio. (1.10) 

If the flow is steady, (1.9) implies that 9 is position independent. Consequently, as we follow a streamline 
through the orifice into the core of the jet, we must have the classical Bernoulli relation 

imw 2 + p, = const. (1-11) 
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Now the pressure in the jet and that of the adjacent ambient fluid must be the same, or the fluid would 
move sideways. This requires the value of \i in the jet to equal /j, Q . Thus the velocity of the fluid in the jet 
is given by 

imv 2 =/ii-^ , (1-12) 

a familiar classical result. We find then that 

-h^j.(6jet ~ ^ambient) = ^mv 2 . (1.13) 

This accumulating phase slip must be accounted for by the passage of vortices in the shear flow bounding 
the jet [11]. 

To verify this we note that the shear flow has vorticity v per unit length. Since 9 winds through 2ir as 
we encircle a vortex, the quantum of circulation is 2nh/m. There are therefore mv/2irh vortices per unit 
length in the jet boundary. Each of these vortices will be convected with the velocity field due to all the 
other vortices, and, by arguments familiar from the calculation of the force on the windings of a solenoid, 
this velocity is v/2. Thus smoke-ring- like vortices are shed at the orifice and convected downstream at a rate 
of ^mv 2 /2nh vortices per second. Each vortex allows a phase-slip of 2ir, so the phase slip between points 
inside the jet and those outside accumulates at a rate \mv 2 /%, in comfortable conformity with (1.13). 

This result is an example of Anderson's general relation, proved in Appendix B of [11], between the 
difference in the Bernoulli constant at two points and rate of passage of vortices across any curve connecting 
them. 

3 Numerical Procedure and Results. 

For the purpose of our numerical work we use the Gross-Pitaevskii equation in the following form. 

id t i> = - v 2 ^ + (M 2 -l)V>- (3.1) 

This corresponds to redefining the units for x and t, and the normalization of ip. In these new units the 
speed of sound c is equal to \[2. As we mentioned earlier, we confine our study to two space dimensions. 

We implemented the time evolution using a stabilized leapfrog algorithm of second order accuracy. We 
used a 200 x 200 lattice with the lattice spacing in spatial directions chosen to be equal to 0.2 and the 
spacing in time direction equal to 0.004. Evolution was stable with these parameters and the total density 
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was conserved to within 1-2 percent during the evolution. The calculations were performed on SPARC 10 
workstation at ITP Santa Barbara and on HP735 workstation at IOP Bhubaneswar. 

We have studied three different cases corresponding to various geometries and sizes of the orifice. Figures 
la, lb and lc show the contour plots for the initial configurations of ip corresponding to these three cases. 
The order-parameter phase, 6, is equal to zero everywhere for these initial configurations. A fixed boundary 
condition ip = is maintained at x = 0. At the other boundary in x direction (x = 40 in our simulations) free 
boundary conditions are used, while periodic boundary conditions are used in the y direction. The midpoint 
of the orifice is at x = 20. Figure la shows the case of orifice with sharp edges where ip — is kept fixed on 
a line along the y axis (at x = 20), leaving the opening for the orifice. Figure lb shows the case when the 
edges of the orifice are semicircles, and Fig. lc shows orifice, again with rounded edges, but with a wider 
opening. Dense regions of contours correspond to the healing length in which ip grows from the value ip = 
to the value \tp\ = 1. 

We begin by using a relaxation method to obtain initial solutions of the G-P equation with the boundary 
conditions described above. We evolve an initially chosen configuration using the diffusive equation 

a t v = v 2 V<-(IVf (3-2) 

After about 2000 time steps tp( x ) converges to a solution of the time independent G-P equation. These 
static solutions are the configurations shown in figures la, lb and lc. They become the initial condition for 
ip for the time evolution using equation (3.1). 

We tried several procedures for generating a superflow. Our first attempt utilised a piston moving along 
the positive x axis from x = 0. The piston was given a small uniform acceleration until it reached a suitable 
velocity. Shedding of vortex-antivortex pairs was observed using this method. Unfortunately, because of the 
finite size lattice, a constant piston velocity led to an increasing average velocity for the superflow through 
the orifice. This method, therefore, was not very suitable in determining the final asymptotic velocity of the 
superfluid through the orifice. 

The second method consisted of the addition of a source term of the form iwip/\ip\ 2 to the G-P equation. 
Here, w is the strength of the source which is taken to be zero initially and increased slowly to a suitable 
maximum value during the course of the simulation. The source was taken to be non-zero along the line x 
= 0.4. This greater control allowed by this procedure provided a constant average superfluid velocity. The 
results reported in the following paragraphs make use of this method. 
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While we used equation (3.1) in most of the region we found it convenient to use the diffusive equation 
(3.2) in the large x region. This helped to damp out sound waves so that the vortex shedding could be 
studied in detail for large time durations. For the orifice with sharp edges, equation (3.1) was used for region 
between x=0 to x=30, while equation (3.2) was used between x=30 to x=40. For orifices with rounded edges, 
these regions were respectively x=0 to x=35 and x=35 to x=40. Since the dissipative equations were used 
for regions sufficiently far from the orifice, in the direction of the flow, it did not interfere with the process 
of vortex shedding. 

Figures 2, 3 and 4 show the details of the process of vortex shedding for the three types of orifices. In 
these plots, solid lines are used to plot contours of p and dashed lines are used for 9 contours. As described 
above, bubbles are pulled off the edges of the orifice. These form the cores of an vortex-antivortex pair 
which detaches from the edge and is swept along the flow. It is interesting to note that, for the orifice with 
rounded edges, vortex shedding does not occur at the point of smallest opening, but further downstream. 
This is not uncommon in boundary layer separation phenomena. It would be interesting to investigate what 
factors determine the point of vortex shedding for such cases. 

Figures 5a, 5b, and 5c show the plots of the average velocity of the fluid flow through the different orifices 
as a function of time. We see that the flow increases uniformly, reaching, on the average, a critical velocity 
asymptotically. This critical velocity is roughly 0.4 (recall c = 1.414) for all three cases, irrespective of the 
shape of the edges of the orifice, or the size of the opening. We see that the fluid velocity first increases, and 
then there is a sudden decrease in the velocity. Looking at the contour plots at these times we see that the 
sharp decrease in the velocity coincides with the shedding of a vortex-antivortex pair. Actually, for rounded 
edges (figures 4b and 4c), the vortex shedding slightly precedes the drop in the fluid velocity. This is not 
unexpected as the vortex is produced downstream from the midpoint of the orifice, while the fluid velocity is 
calculated at the midpoint. It will take a while before the vortex shedding can affect the midpoint velocity. 
This further supports the idea that the drop in fluid velocity is caused by vortex shedding. Each successive 
large jump in the fluid velocity corresponds to a vortex-antivortex pair being shed. 

An intuitive understanding of this sudden jump in the fluid velocity can be given in the following way. 
As the core of the vortex buldges out of the edge, it effectively narrows the orifice. This leads to an increase 
in the fluid velocity through the orifice which in turn further lowers the pressure. Eventually the vortex is 
shed and the backwards fluid motion around the vortex tends to reduce the net flow through the orifice. A 
significant part of the kinetic energy of the flow is transferred to the vortex. 
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4. Discussion 

The vortex shedding process exhibited in this paper is probably not directly relevant for interpreting 
the experiments. The critical velocities reported in refs. [4,5] seem too low to cause cavitation. It is more 
likely that some form of quantum tunneling pre-empts the classical mechanism. It is, however, well worth 
obtaining a better understanding of the classical mechanism. The spontaneous shedding of vortices seen 
in our computation demonstrates that, at least for G-P fluids in two dimensions, the tunneling barrier to 
creating a vortex disappears above some critical velocity. This critical velocity needs to be estimated if the 
purely classical effect is to be distinguishable from the more interesting macroscopic quantum tunneling. 

How does the vortex tunnel? Homogeneous quantum nucleation [12,13] is not possible when there is 
no normal fluid to select a preferred reference frame, but there are several competing channels that use 
the boundary to break galilean invariance. One possibility is to create a complete vortex ring surrounding 
the orifice. This, however, requires a large action. A more popular scenario is the quantum nucleation of 
a semicircular segment of vortex terminating on the boundary [14,15,16,17]. This segment is then swept 
into the flow where it grows by extracting energy from the bulk motion. Because this likely scenario breaks 
the axial symmetry, two-dimensional computations are not adequate for determining the classical critical 
velocity — if it exists. (Both the thin vortex models of [14,15] and the static G-P analysis in [16] suggest 
that a tunneling barrier persists even for large flow velocities.) It is therefore essential to perform a full 
three-dimensional time-dependent computation for an orifice with a suitable asperity on which the vortex 
can form. 
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Figure Captions. 

Fig 1. Initial configurations. The solid lines are contours of density. 

Fig 2. Early stage of flow. The solid lines are contours of density. Dashed lines are contours of order-parameter 
phase. The heavy dashed lines are branch cuts where the order-parameter phase jumps through 2tt. 

Fig 3. Same as Fig 2., but at later time. 

Fig 4. Same as Fig 3., but at later time. 

Fig 5. Flow at midpoint of aperture. Each sharp dip in the velocity corresponds to a vortex-antivortex pair 
being shed. 
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Contour plot at t = 0.0 




Fig. 1b 



Contour plot at t = 79.0 




Contour plot at t = 68.0 




Fig. 2a 



Contour plot at t = 82.0 




Contour plot at t = 72.0 




Fig. 3a 



Contour plot at t = 84.0 




Fig. 4b 



Contour plot at t = 0.0 
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Contour plot at t = 58.0 
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Contour plot at t = 62.0 




Contour plot at t = 64.0 
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Contour plot at t = 76.0 




Fig. 4a 



